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Abstract 

We study the anomalous density in an ultra-cold trapped bose gas in a variational 
framework, for both zero and finite temperature. We show that it is finite in ID, 
while it is logarithmically and linearly divergent in 2D and 3D. The renormalization 
that we adopt is more reliable and compatible with the variational scheme. The 
main outcome is that the anomalous and non condensate densities are of the same 
order of magnitude and should therefore be treated on an equal footing. 
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1 Introduction 



In a previous paper [1], using the time dependent variational principle of Balian and 
Veneroni[2], we have derived variationally a set of dynamical coupled equations (which 
we called TDHFB for Time-Dependent Hartree-Fock-Bogohubov) for dilute trapped bose 
gases below the transition. These equations were shown to generalize in a consistent way 
the Gross-Pitaevskii equations[3] in that they introduce a non trivial coupling between 
the order parameter, the thermal cloud and the anomalous density. 

In order to apprehend better this coupling and to understand the role of the anoma- 
lous density, we extended the formalism to the non-local case [4, 5], thus obtaining an 
intrinsic dynamics of the thermal cloud and the anomalous average which has never been 
written down before (an exception is perhaps the paper of Chernyak et al. [6] which 
discusses a somewhat similar set of equations using the generalized coherent state repre- 
sentation). We focused in these last two papers on the anomalous density owing to its 
importance to account for many-body effects. Nonetheless, the subsequent analysis was 
mainly performed in the low momentum regime, therefore leading to an effective theory 
which takes into account only the low lying modes. In this respect, the divergences of 
the anomalous density or the many-body T-matrix which are known to occur at high 
momenta were hidden in our treatment. 

Despite the fact that the resulting equations were greatly simplified in the static 
case, leading to useful conclusions about the anomalous density, it was not obvious that 
the inclusion of higher modes will maintain these results, especially at finite temperature. 
In particular, it was observed [7] that the coupled dynamics of the atoms in and out of the 
condensate induces an effect which is highly needed in order to go beyond the HFB-Popov 
approximations. Indeed, the discrepancies between the HFB-Popov computed properties 
and the experimental resuts that show up at temperatures as high as 60% of the transition 
temperature seem to be related to the absence of a feedback effect of the dynamics of both 
the thermal cloud and the anomalous average on the condensate density [8]. 

In this paper, we would like to go beyond the approximations discussed in [4, 5] 
by explicitely taking into account the high lying modes. This will obviously lead to 
self-consistent equations exhibiting UV divergences that have to be regularized. 

In this context, we note that in reference [6], the diverging terms misteriously dis- 
appear from the final equations. The subsequent analysis performed there is thus of a 
very limited validity. In [9] however, the A potential technique [10] is developed in order 
to obtain a consistent HFB theory. This formahsm is however not fully variational and 
endows many intricacies such as the parameter A itself. In our case, we will use instead 
the Wigner representation[ll, 12, 13] which is best suited for nonuniform systems and 
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which will allow us to handle the divergences in a variational framework. We will see that 
we arrive at the same renormalization scheme as [9] in a more reliable way. Moreover, we 
obtain important quantitative estimates for the anomalous average and the non conden- 
sate density leading to the fact that they are essentially of the same order of magnitude. 
We therefore arrive at the conclusion that what is called in the literature the HFB-Popov 
approximation [7, 8, 14, 15] is a quite hazardous and highly inconsistent approximation [16]. 

The paper is organized as follows. In section 2, we recall the main steps in the vari- 
ational derivation of the time-dependent Hartree-Fock-Bogoliubov equations for a system 
of trapped bosons with general two-body interaction and study some of their properties. 
Then we speciahze to the contact potential case which allows us to show that the HFB- 
Bogoliubov-De Gennes (HFB-BdG) equations are just the quasi-homogeneous hmit of our 
general equations. In section 3, we discuss the zero temperature static solutions. We focus 
in particular on the anomalous density owing to its importance to account for many-body 
effects. We derive analytic expressions for both the anomalous and the noncondensate 
densities for one, two and three dimensional systems. The anomalous density is shown 
to be finite for D = 1 while it is logarithmically and linearly divergent for D = 2 and 3. 
Our renormalization scheme is not only more reliable and consistent with the variational 
approach, but it leads to the same condition of renormalization found with the pseudo- 
potential methods. Most importantly, the anomalous and the non condensate densities, 
although small at zero temperature, are found to be of the same order of magnitude, 
leading therefore to the conclusion that it would be very hazardous to neglect one of 
them while keeping the other such as what is performed in the so (wrongly) called Popov 
approximation. 

In section 4, we evaluate the leading terms in the low temperature approximation 
and observe that our conclusions are not altered by the finite temperature corrections. 
Finaly, we present some concluding remarks and perspectives. 

2 The TDHFB Equations 

Let us first recall how the generalized TDHFB equations may be derived in a variational 
framework. The Balian-Vcncroni variational principle requires the choice of a trial space 
for the density operator. We will consider a gaussian time-dependent density-like operator. 
This Ansatz, which belongs to the class of the generalized coherent states, allows us 
to perform the calculation (since there exists a Wick's theorem) while retaining some 
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fundamental aspects such as the pairing between atoms. Let us first define the quantities 

$(r) 



m r, r 



n*(r',r) = {-ip+{r)-ip{r')) - $*(r)$(r'), 
m(r',r) = {'ijj{r)'ijj{r')) - $(r)$(r'), 



(2.1) 



where il^{r) and (r) are the boson destruction and creation field operators satisfying 
the usual commutation rules [ip(r),^~^(r')] = 5{r — r'), [ip(r),^(r')] = [^~^ (r) , {v )] = 0. 
These quantities completely characterize the gaussian density operator and, in our lan- 
gage, are akin to an order parameter and to one-body correlation functions. The local limit 
of the latter are just the noncondensate density limj./^j. n(r, r ) = n(r) and the anomalous 
density lim^'^^ m(r, r ) = m(r). Let us consider a second quantized hamiltonian H of the 
form: 

H = j^^+(r)i/^P(r)^(r) + ^ ^+(r)^+(r>(r, r>(r>(r), (2.2) 
where V{r, r) is the interaction potential and if^P(r) the single particle hamiltonian 



^^'(r) = -^Ar + Kxt(r)-/x. 
Zm 



(2.3) 



T4xt(r) and /i are respectively the external trapping field and the chemical potential. 
The stationarity conditions of the Bahan-Veneroni action-like functional lead to a set of 
dynamical equations which couple the various densities: 



ihn{r, r') 



+ V{r, r') |$(r')p<l>(r) + $*(r')m(r, r') + $(r')n(r, r') + $(r)n(r', r' 



n[r, r 



+ Jj." V{y , r ) a{r , r )n(r, r ) + a(r , r )n(r, r ) + o(r , r )m (r , 

„ T7-/ "\ [ / II \ ^ / II l\ /II "\ — / '\ T=k/ II \ ^ f II l\ 

— Jj." V{r,r ) a(r, r )n(r ,r)-|-a(r ,r )n(r, r ) + o (r, r )m(r ,r) 
= i/^P(r) + H^P{r')] m(r, r') 

+ Jr" ^ ) •> ^ r ) + a{v , r )m(r, r j + o(r . r ){7i (r. r ) + o[y — v )) 

„ , " \ \ / " \ ~ / ' " \ / " " \ ~ / ' \ 7/ II \ ~ / II I 

+ Jj." K (r, r ) a(r , r)m(r , r ) + a(r , r )m(r, r j + 6(r, r )n(r , r 



(2.4) 



where 



and 



a(r, r') = a*(r', r) = (V''^(r)V'(r')) = fi(r, r') + $*(r)$(r') 



6(r, r ) = 6(r , r) = (■?/'(r)-?/'(r )) = m(r, r ) + $(r)$(r ). 
In the case of a contact potential V{y,y') = gS{r — r ), the equations (2.4) read 

ih^r) = [//'P(r) + g(nc(r) + 2h{r))] $(r) + gm(r)^*(r), 



(iJ'^P(r) + 2gn{r)) - (iJ^P(r') + 2gn{r' 



+9 



b{r\ r')m*(r, r') — 6*(r, r)m(r, r') 



(2.5) 
(2.6) 
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ihrh{r, r') = [{H^^r) + 2gn{r)) + (h^^t') + 2gn{r'))] m(r, r') 
+9 [b{^' ,r'){5{r - r') + n*{r, r')) + b{r, r)n{r, r')] , 

where nc(r) = |$(r)p is the condensate density and n(r) = ndr) + n{r) the total density. 
For convenience, we have omitted the time dependence. We notice that even for a contact 
potential, the equations (2.6) and (2.7) are still complicated. Indeed, the presence of sec- 
ond order derivatives (in iJ^P) and the 6 term in (2.7) do not allow for simple local hmits. 
This is precisely the term which misteriously disappear in the derivation of Chernyak et 
al.[6\ using the coherent state representation. We will see in section 3 that this term is 
responsible for the divergences while its structure allows for a correct renormalization of 
the coupling constant. 

On the other hand, this set of equations was also derived by [14, 15, 17] in a local 
form by a quite different approach. In our case however, these equations can be related to 
the well-known HFB-BdG equations only in the quasi-homogeneous case. Indeed, upon 
setting 

^(r, r') = -1 Efc(l + 2n,) {Uk{r)Vk{v') + C/fe(r')14(r)) , 
ri(r, r') = {n,U*,{v)Uk{v') + (1 + n,)V,*{v)V,{v)) , 

where Uk and are linearly independent space functions and rik is the occupation prob- 
ability of the mode k which is given at equihbrium by the Bose-Einstein distribution 
function, we readily get from (2.6) and (2.7) 

ihUkir)= [H^P{r) + 2gn{r)]Uk{T)-gb{T,T)V:ir), 
ihVkir)^ [H^P{r) + 2gn{r)]Vk{r)-gb{r,r)U^{r), 

which, together with (2.5) form the so-called (generahzed) time dependent HFB-BdG 
equations [7, 11, 14]. Hence, our equations can be considered as a variational extension 
of the time dependent HFB-BdG equations for non-homogeneous situations. 

Among the interesting properties of the TDHFB equations, wc can quote the con- 
servation of the energy and of the total particle number. Moreover, the "Heisenberg" 
parameter, which is defined as 

7(r,r') = ^^^ ((5(r-r") + 2fi(r,r")) (5(r" - r') + 2fi(r", r')) - 4 m*(r, r")m(r", r'), 

(2.10) 

and which can be written in the more readable compact form 

/(r,r') = (5(r-r') + 2n(r,r'))^-4|m(r,r')|^ (2.11) 

is also a constant of motion at zero temperature. Indeed, for T = 0, /(r, r ) = 6{r — r') and 
/j." n*(r, r ")rn(r ", r') = /j." m(r, r ")n(r ", r ). One can easily check that these properties 
remain true during the evolution if they are initially fulfilled. This conservation law is not 
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however new since it is intimately related to the unitary evolution of the single particle 
density matrix [1, 2, 18]. 

The equation (2.7) involves an explicitly diverging term (for r = r') which should 
be regularized. To this end, one may absorb this diverging expression in a redefinition 
of the thermal average[7, 8, 19] or, as performed by [20, 21], one may use the pseudo- 
potential method [22] to renormalize the coupling constant. A more rigorous approach is 
the A-potential method first discussed in[9, 10]. Neverttheless, these approaches suffer 
from conceptual difficulties in that they are based on substr action schemes, which are not 
compatible with the variational technique. Another way of proceeding, which seems to us 
much more interesting, since it allows to remain in the variational framework, is to resort 
to a better representation which would "wash out" these singular terms. One of the best 
candidates is the Wigner representation[ll, 12, 13]. To this end, let us define the Wigner 
transform for any function A{r,r') as 

A(r,r') = I^AwC-^,k)exp[-ik.{v-r)], (2.12) 

and its inverse 

^^'(R', = I ^(R- + - r/2) exp (ik.r), (2.13) 

where = ,J'\ti I ^^k. In the hmit r = r', the equations (2.6) and (2.7) take the form 
[Ztt) 



.t2 

ihfiwi^,^) = -77^ k.Vrnw/(r, k) 



+g ($2(r) + 'm{r))'m^{r, k) - ($*^(r) + m*{r))mw{r, k) 



(2.14) 



ihrhwir, k) = 2 



H'P + 2gn{r) + 



^ miy(r,k)+<7($2(r)+m(r))(l+2niy(r,k)). (2.15) 
2m 



In the following sections, we will provide the static solutions to these equations 
both at zero and finite temperature. This will not only yield quantitative estimates of the 
anomalous effects, but also lead to some important remarks concerning the behaviour of 
the anomalous density with respect to the dimensionality of the problem. 



3 Anomalous Effects at T = 

In the static case, and for a neutral gas, one may consider without loss of generality that 
m and n are real space functions. The equation (2.14) translates into 

k.VrnvK(r,k) = 0, (3.1) 
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which means in particular that there is no current associated with the non condensate and 
that the condensate current is therefore conservative. On the other hand, the equations 
(2.5) and (2.15) yield 

0= [i/'P + ^(nc + 2fi + m)]$(r), (3.2) 

= 2 ( i/^P + 2gn{r) + ^ j rhwir, k) + g{n, + m)(l + 2hw{r, k)). (3.3) 

Next, one may use the gradient expansion method[ll, 12, 13] to show that riw and rhw 
write simply 

2mw{r,k) = -(1 + 2nA;) sinh f7(r, k), 
l + 2fiH^(r, k) — (1 + 2nfe) coshcr(r, k), 

where we recall that rik is the Bose-Einstein distribution function at equilibrium. The 
Wigner transform of the Heisenberg parameter Iw (2.11) takes the familiar form 

= l + 2nfe = coth(/3efc/2), (3.5) 

(3 being the inverse temperatures. The quantities are the s.p. energies obtained by 
solving the stationary Gross-Pitaevskii equation (3.2): 

efe = — \-Ve^t - 11 + g{n + n + m). (3.6) 

zm 

Corrections to the expressions (3.4-3.5) will measure how far the gas is from homogeneity [11] . 

The spectrum (3.6) is manifestly gapless. Indeed, in the homogeneous long wave- 
length limit, the Hugenholtz-Pines theorem[23] is automatically satisfied (see the end of 
this section). In contrast, in the standard HFB framework, one has to slightly modify the 
equations in order to obtain a gapless spectrum (see e.g. [7, 8, 14, 24]). Similarly, in the 
A-potential approach [9, 10], a zero gap condition is imposed in order to render the theory 
consistent and this leads to the correct equation of state. 

In order to analyze the anomalous effects, let us derive workable expressions for n 
and m. To this end, we follow the method depicted in [13]. The key point shown there 
is that, up to order fi^, the Wigner transform of the product of two s.p. operators is 
the product of their Wigner transforms [11, 12]. The equation (3.2) may thus be written 
formally in a classical form 

i^i-H,i)^ = 0, (3.7) 
with the self-consistent hamiltonian Hd = /2m + V^cxt + g{i^ + n + fh). The local mo- 



mentum p must satisfy the on-shell condition p = po(r) = •\/2m(/x — V^xt ~ fl'(^ + n + m)). 
Using (3.4), we get the intermediate result 

tanh a r, k = — ^ 3.8) 

+ Kxt - A* + 2yn + 



2m 2m 
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which yields the desired expressions 



m 



for n and m: 

gijic + fh) 



(27r) 



^ V 2m 



2l2 



, 2m 



= coth(/3efe/2), (3.9) 



— 2gm 



n 



(2^ 



2 7,2 



■1 + 



2m 



+ g{nc - fh) 



\ 



2l2 



'h'k 



2l2 



: coth (/3efe/2) 



— 2gfm 



(3.10) 



2m J \ 2m 

What is remarkable on these expressions is that the denominator is clearly not of the 
Bogoliubov type. This is the price to be paid in order to maintain not only the consistency 
of our variational approach but also the gapless nature of the spectrum. To recover a 
Bogoliubov-type spectrum, one may set m = in the integrand, which is reminiscent to 
the so-called Popov approximation. But this will be shown to be highly inconsistent in a 
variational framework. 

It is also important at this level to notice that the on-shell condition is somewhat 
equivalent to the local density approximation in the sense that the effect of the confining 
potential is only present in the temperature dependent term. 

Furthermore, since the previous expressions can be generalized to any dimension D 

by simply making the replacement / " . ■> — )■ / r, , one can see that they are both 

(27r) (27r) 

free from IR divergences. On the other hand, while n is finite at high momenta, we can 
clearly see that m is UV divergent for D = 2 and D = 3. It is however finite in the 
one- dimensional case. 

In the latter case and at zero temperature, a typical length scale (the interaction 
length) is ai — /mg. A straightforward calculation then yields: 



n/ric 



m\ n, 



^arctanh 



1 - 



1 



ml Ur 



\m\/nc = ^y'l - \m\/nc (log (1 + ^Jl- jm/nd) - log y^jm/nd^ 



(3.11) 



(3.12) 



where we have introduced the dimensionless gas parameter 7 = 1/ ^/aln^ which measures 
the ratio between the interparticle length l/ric and the interaction length ai. These 
simplified expressions are obtained under the assumption that \rh\/nc << 1. In the 
general case, they will yield the elliptic functions K and F. Using this last assumption, 
we may further simplify the equations (3.11) and (3.12) to get: 



n/ric ~ ^ log {l/yj\rh\/nc), 



(3.13) 
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|m|/ne~^log(l/^|m|/nc). (3.14) 

We note on these expressions that the anomalous fraction is half the quantum depletion. 
Although they are both small, it would be hazardous even at zero temperature to neglect 
the anomalous average while maintaining the non condensate density. Moreover, the 
equation (3.14) does not admit solutions unless 7 << 1 which is the the weak coupling 
condition aiUc » 1 [25]. In this case, the explicit solutions may be written as: 

^=^^^^v^(l-V7)- (3.15) 

However, the next to leading correction of order 7 is not relevant since it seems to be 
related to three-body losses which are not considered here[26] . Hence, up to order ^/7, h 
and m are of the same order and scale as y^. 

In order to extend the calculation to more than ID, one must handle carefully the 
UV logarithmic (2D) and hnear (3D) divergences that appear in the integral (3.9). 

Let us begin by the 2D case. First of all, the expression (3.10) can easily be evaluated 
to yield 

n=——^nc[l^ , (3.16) 



where q = ^ is very close to 1 at low temperatures. On the other hand, the integral 
(3.9) must be regularized. Upon introducing a cutoff K, we obtain for high momenta 

^)log(2X/Xo) 
^ ^-¥Ur-^ ^c, (3.17) 



where Kq = ^^2^ rig. We may now renormalizc the coupling constant by setting gji = 
g log {2K/ Kq) and introduce the dimensionless gas parameter 7 = mgR/2TTh^ to obtain a 
finite expression for m: 

in = ^ — ric- (3.18) 

1 + 7 

It should be emphasized that this renomalization of the coupling constant is analogous to 
what is obtained in [9] (see equation (50) of this reference). Furthermore, our cutoff K is 
just the parameter A introduced there. The final expressions for m and n (replacing in 
the latter g by gu and setting ? ~ 1) read 



- " ^ (3.19) 



Uc Uc 1 + 7 

which show that h and m are also of the same order but now scale as 7 
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In 3D, although the calculations are mainly the same, their outcomes are quite 
different. Indeed, we obtain for h and rh the expressions 

n^-^[a{n, + m)f\ (3.20) 

fh = —aK{nc + m), (3-21) 

where X is a high momentum cutoff and a is the "bare" s-wave scattering length, which 
satisfies aK — —ih/{nc + rh). This is the right renormalization condition since upon 
redefining the coupling constant by 

we get Qr = g{\. + ^) which redefines both the scattering length a — >■ = mgii/ATrh'^ 
and the noncondensate density 

n = —= [aRUc] . (3.23) 

The previous renormalization condition was also derived in [9, 10] based on the A-potential 
and by [7, 8] using a somewhat ad-hoc prescription. In this sense, the approaches devel- 
oped in the previous references are compatible although they seem quite different. It is 
however remarkable that (3.23) can be directly deduced from (3.20)by setting m = 0, 
which is also compatible with (3.21) if we set a — J2s=ils/ . In this case, we get the 
finite part of fh 

rh = ^ne. (3.24) 

1 + 71 

Hence, unlike the ID and 2D cases, rh and h are not of the same order of magnitude. The 
former takes values compatible with 0, as in the HFB-Popov formalism. 

Finally, one may derive the equation of state in each case using the Hugenholtz-Pines 
theorem which relates the chemical potential to the diagonal and off-diagonal self-energies 
^Eii and ^Ei2 hy /i — HT^n — hl^i2- The latter are given by 

= ^ (3.25) 

where To is the D-dimensional T-matrix. In the case of a contact potential, is just 
the bare coupling constant g. 
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4 Anomalous Effects at Finite Temperature 



In order to compute the finite temperature corrections, one may rewrite (3.10) in the 
form: 



n = uq 



2 1 (2.)^ 1/ 2 2 X / 2 2 (coth(/3e./2)-l). (4.1) 

+ 2gnc — 2gm 



^ y 2m '^J \ 2m 

where, for low T, the first term may be approximated by its value at T = 0, hence the 
subscript 0. The second term may be computed as follows. At very low temperatures, 
the integral is dominated by the behavior of its integrand for k = ko = po/fi, Po being the 
local momentum given by (3.7). Upon setting 6n = n — no, we obtain: 

c r(n/o\ \- Q(nco — "mo) 

Sn ~ ' / ^ 2m \-D g-/9(Vext-M+ff(no+no+mo))^ ^4 2) 

where S*! = 1, 5*2 = 27i, S3 = in, T is the Euler function and A is the thermal De 
Broglie wavelength. The subscript in the densities means that they are evaluated at 
zero temperature. Proceeding in the same way for m, we get 

The previous expressions show that the finite temperature corrections although being 
exponentially small are of the same order for both m and n. Our conclusions concerning 
the order of magnitude of these two quantities remain therefore valid at finite temperature. 

One may further simplify these expressions owing to the equation of state and to 
the results obtained in section 3. The main point here is that the temperature dependence 
of the corrections, beside the exponential factor, is T^/^. 



5 Conclusions 

In this paper, we have been mainly interested in obtaining quantitative estimates of the 
so-called anomalous average for bose condensed gases. Indeed, owing to its importance to 
account for many-body effects, it was not clear for us what is the status of this quantity 
especially in the HFB formalism. 
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To this end, we preferred to use from the beginning a consistent formahsm. The 
time-dependent variational principle of Balian and Veneroni appeared as the best candi- 
date. Indeed, not only does it lead to approximate dynamical equations by a judicious 
choice of the trial spaces, but it maintains also the consistency of the approach along the 
whole calculational process. In this context, we obtained first a set of time-dependent 
equations for the order parameter and the one-body correlation functions, which were 
shown to generalize the HFB equations in the non-local, non-homogeneous case. More- 
over, using the Wigner representation, we were able to take the correct local limit, showing 
that the anomalous density and the non condensate density are the local limits of the pre- 
vious correlation functions. This demonstrates in particular that the HFB-BdG equations 
are only valid in the quasi-homogenous case. Hence, our formalism allows to analyze the 
anomalous effects in more general situations. 

It is to be recalled that we have previously performed a similar analysis in the low 
momentum regime and obtained quite interesting results about the anomalous density. 
The question that naturally arises was whether these results remain true when we take 
into account the high lying modes. This question becomes crucial at finite temperature. 

Our analysis shows that the anomalous density and the non condensate density 
are of the same order in one and two space dimensions at zero and finite temperature. 
Therefore, neglecting only the former is a quite hazardous approximation. In 3D however, 
the situation is much more subtle and requires a careful analysis. A preliminary result is 
that the anomalous density is effectively small compared to the non condensate density. 

Moreover, using a renormalization scheme compatible with our variational approach, 
we obtained in a more reliable and rigorous way, finite expressions for the anomalous 
density and the non condensate density similar to what is obtained in the literature. 

Many unanswered questions remain to be examined. One of the most important 
is whether the previous results depend on the nature of the interactions. Indeed, the 
contact interaction that we adopted in this paper is the simplest and care should be taken 
especially for long range interactions. We plan to publish in the near future a similar 
analysis for dipolar gases. 

Moreover, it has been shown elsewhere that the HFB approximation consist in fact of 
a resummation of ladder diagrams[9] which opens the door to the well known speculations 
about the relationship between the variational and the perturbative approaches [2 7]. 

The time dependent renormalization is also of great importance. In our variational 
formalism, we feel that this can be done more easily since we have shown elsewhere [28] 
that the time dependent mean field equations are renormalizable once they have been 
renormalized initially. 
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